Preguntas de investigación

  1. Describir autocorrelación espacial, autocorrelación espacial positiva, y autocorrelación espcial negativa
    • Autocorrelación espacial: la autocorrelación espacial nos ayuda a identificar el nivel de de correlación que tiene el valor de una variable consigo misma en el espacio, en palabras más simple ayuda la similitud de objetos u observación en un plano espacial (R Spatial, 2019)
    • Autocorrelación espacial positiva: es cuando existe un cluster u agrupamiento de valores similares en el espacio, es decir cuando los valores vecinos son similares, este es el tipo de autocorrelación más común (GISGeography, 2014).
    • Autocorrelación espacial negativa: contrario a la autocorrelación espacial positiva, la negativa se da cuando valores disimilares se agrupan en el espacio, es decir los vecinos tienen caracteríssticas diferentes (GISGeography, 2014)
  2. Describir los conceptos de autocorrelación espacial y no estacionariedad en un contexto de análisis espacial de datos
    • Autocorrelación espacial busca mostrar el nivel de similitud de un objeto con sus vecinos en el espacio, la no estacionaridad espacial se da cuando la variable de interés no es constante en la región estudiada, también conocido como heterogeneidad espacial, cuando esto sucede existe una violación al supuesto de estacionaridad, que afirma una constante en las propiedades estadísticas de una variable en el espacio (GIS&T, 2016)
  3. Describir al menos 3-5 diferencias entre análisis exploratorio de datos (EDA) y análisis exploratorio espacial de datos (ESDA)
    • EDA no investiga explícitamente el componente espacial de un set de datos, mientras que ESDA sí lo hace (Abdishakur, 2019)
    • Al momento de hacer la visualización de datos en el EDA no se hace un discriminación de locación geográfica mientras en el ESDA sí, por lo cuál mucha de su representación se hace directamente en mapas (Abdishakur, 2019)
    • Los análisis ESDA suelen aplicarse en ramas como la geografía, la ecología y el crimen, por su alta necesidad de tomar en cosideración el espacio como una variable de alto peso, los análisis EDA son más usados en ramas como los negocios y ciencias sociales donde la investigación se centra más en insights de patrones no espaciales (Dall’erba, 2009).
    • En cuanto a la estructura de los datos para el EDA no es necesario el componente espacial, mientras que para ESDA es necesario contar con datos georeferenciados y no simplemente el nombre de la localidad (Abdishakur, 2019)
  4. Escribir al menos 3-5 diferencias entre GWR y GRF
    • La principal diferencia entre un GWR y GRF es que debido a la naturaleza de boostrapping de un GFR es más difícil que este modelo sufra de overfitting pues relaja los supuestos de las estadísticas gaussianas tradicionales (R: Geographically Weighted Random Forest, 2019)
    • El GWR se basa en una regresión, mientras GRF se basa en un modelo de machine learning (R: Geographically Weighted Random Forest, 2019).
    • El GWR requiere de una base de datos georeferenciada mientras mientras el GRF al ser un modelo de machine learning puede trabajar con cualquier tipo de variable predictora (Bandwidth selection for basic GWR, 2022).
  5. Describir 3-5 recomendaciones para reducir o eliminar la presencia de autocorrelación espacial en los residuales de un modelo de regresión estimado.
    • Transformar los datos usando el logaritmo o la raíz cuadrada de la variable dependiente, pues ayuda a reducir el peso de los outliers que a su vez afectan la autocorrelación espacial en los residuos (Guélat, 2013).
    • Así como se puede transformar los datos usando log o raíz cuadrada se pueden directamente quitar los outliers o modificarlos con la media (Guélat, 2013).
    • Reevlauar la elección de variables de entrada hasta que ya no exista autocorrelación significativa (ArcGIS, 2019).

Análisis Exploratorio de Datos (EDA)

data(columbus) ### dataset
columbus_shp <- readShapePoly(system.file("etc/shapes/columbus.shp", package="spdep"))
## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among others
### shapefile
#col.gal.nb <- read.gal(system.file("etc/weights/columbus.gal", package="spdep"))
### matriz de conectividad espacial pero requiere calcular los pesos espaciales.
##eliminar columnas con nombres dupliccados
columbus <- columbus[, -c(1, 2)]

Descripción de variables de interés

columbus <- columbus %>% select(HOVAL, INC, DISCBD, CRIME, PLUMB)

Se seleccionarone estas variables según la descripción mostrada, pues de toda la base de datos fueron las características consideradas más relevantes para la variable dependiente que en este caso será HOVAL

HOVAL: valor de la vivienda (en $1,000)

INC: ingresos del hogar (en $1,000)

CRIME: robos residenciales y robos de vehículos por cada mil hogares en el barrio

PLUMB: porcentaje de unidades de vivienda sin plomería

DISCBD: distancia a CDB

summary(columbus)
##      HOVAL            INC             DISCBD          CRIME        
##  Min.   :17.90   Min.   : 4.477   Min.   :0.370   Min.   : 0.1783  
##  1st Qu.:25.70   1st Qu.: 9.963   1st Qu.:1.700   1st Qu.:20.0485  
##  Median :33.50   Median :13.380   Median :2.670   Median :34.0008  
##  Mean   :38.44   Mean   :14.375   Mean   :2.852   Mean   :35.1288  
##  3rd Qu.:43.30   3rd Qu.:18.324   3rd Qu.:3.890   3rd Qu.:48.5855  
##  Max.   :96.40   Max.   :31.070   Max.   :5.570   Max.   :68.8920  
##      PLUMB        
##  Min.   : 0.1327  
##  1st Qu.: 0.3323  
##  Median : 1.0239  
##  Mean   : 2.3639  
##  3rd Qu.: 2.5343  
##  Max.   :18.8111

Se pueden observar diferentes estadísticos que muestran la tendencia de las variables, más adelante se confirmará en los gráficos como variables como PLUM no siguen una distribución normal.

plot_normality(columbus, HOVAL, INC, DISCBD, CRIME, PLUMB)

Observando los histogramas de cada variable se pueden ver como varias no siguen una distribución normal, para poder reducir su sesgo se harán transformaciones logaritmicas o de raíz cuadrada a conveniencia de cada varible en particular.

boxplot(columbus$HOVAL)

boxplot(columbus$INC)

boxplot(columbus$DISCBD)

boxplot(columbus$CRIME)

boxplot(columbus$PLUM)

Con los boxplot hechos la presencia de outliers puede confirmar la falta de distribución normal en las variables.

Transformación de variables

columbus$HOVAL <- log(columbus$HOVAL)
columbus$INC <- log(columbus$INC)
columbus$PLUM <- log(columbus$PLUM)

Se realizó una transformación logaritmica a aquellas variables que contaban con falta de normalidad, siendo un cambio directo a la base de datos ya no será necesario especificar su transformación en los modelos.

Análisis Exploratorio Espacial de los Datos (ESDA)

swm_queen <- poly2nb(columbus_shp, queen = TRUE)
rswm_queen <- nb2listw(swm_queen, style = "W", zero.policy = TRUE)
coords <- coordinates(columbus_shp)

Autocorrelación espacial global

moran.test(columbus_shp$HOVAL, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$HOVAL  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 2.2071, p-value = 0.01365
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.180093114      -0.020833333       0.008287783
moran.test(columbus_shp$INC, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$INC  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 4.7645, p-value = 9.467e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.415628778      -0.020833333       0.008391926
moran.test(columbus_shp$CRIME, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$CRIME  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 5.5894, p-value = 1.139e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.500188557      -0.020833333       0.008689289
moran.test(columbus_shp$PLUMB, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$PLUMB  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 5.2488, p-value = 7.656e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.416170997      -0.020833333       0.006931979
moran.test(columbus_shp$DISCBD, listw = rswm_queen, zero.policy = TRUE, na.action = na.omit)
## 
##  Moran I test under randomisation
## 
## data:  columbus_shp$DISCBD  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 8.7904, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.800562617      -0.020833333       0.008731396

Este test es usado para mostrar la autocorrelación global de las diferentes variables de interés.

HOVAL: sí existe presencia de auotcorrelación espacial significativa, débil y positiva.

INC: sí existe presencia de auotcorrelación espacial significativa, moderada y positiva.

CRIME: sí existe presencia de auotcorrelación espacial significativa, alta y positiva.

PLUMB: sí existe presencia de auotcorrelación espacial significativa, moderada y positiva.

DISCBD: sí existe presencia de auotcorrelación espacial significativa, alta y positiva.

La variable de HOVAL es la variable con autocorrelación más débil y la de DISCBD con autocorrelación más fuerte.

Autocorrelación espacial local

columbus_shp$sp_HOVAL<-lag.listw(rswm_queen,columbus_shp$HOVAL,zero.policy=TRUE)
columbus_shp$sp_INC<-lag.listw(rswm_queen,columbus_shp$INC,zero.policy=TRUE)
columbus_shp$sp_CRIME<-lag.listw(rswm_queen,columbus_shp$CRIME,zero.policy=TRUE)
columbus_shp$sp_PLUMB<-lag.listw(rswm_queen,columbus_shp$PLUMB,zero.policy=TRUE)
columbus_shp$sp_DISCBD<-lag.listw(rswm_queen,columbus_shp$DISCBD,zero.policy=TRUE)

En los siguientes mapas se hará la demostración de las variables antes de ser transformadas y después, a fin de mostrar las concentraciones o clústers.

qtm(columbus_shp, "HOVAL") 

qtm(columbus_shp, "sp_HOVAL")

qtm(columbus_shp, "INC") 

qtm(columbus_shp, "sp_INC")

qtm(columbus_shp, "CRIME") 

qtm(columbus_shp, "sp_HOVAL")

qtm(columbus_shp, "PLUMB") 

qtm(columbus_shp, "sp_HOVAL")

qtm(columbus_shp, "DISCBD") 

qtm(columbus_shp, "sp_DISCBD")

Se puede observar como de manera local para todas las variables existen más agrupaciones en el lado este del mapa. También es notorio como aumentan las concentraciones una vez se hace la trasnformación espacial de los datos. sp_DISCBD es la única variable que muestra también fuertes agrupaciones en la parte oeste del mapa.

Estimación de Modelos de Predicción

Modelo no espacial

non_spatial_model = lm(HOVAL ~ INC + DISCBD + CRIME + PLUMB, data = columbus) 
summary(non_spatial_model)
## 
## Call:
## lm(formula = HOVAL ~ INC + DISCBD + CRIME + PLUMB, data = columbus)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38117 -0.23022 -0.06776  0.25153  0.91473 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.23809    0.50544   6.406 8.52e-08 ***
## INC          0.11262    0.15501   0.726  0.47139    
## DISCBD       0.11446    0.05183   2.208  0.03247 *  
## CRIME       -0.01175    0.00438  -2.683  0.01025 *  
## PLUMB        0.04600    0.01418   3.243  0.00226 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3078 on 44 degrees of freedom
## Multiple R-squared:  0.5356, Adjusted R-squared:  0.4934 
## F-statistic: 12.69 on 4 and 44 DF,  p-value: 6.008e-07
AIC(non_spatial_model)
## [1] 30.31066

Muestra ser un modelo significativo con todas las variables siendo de igual manera significativas, menos la de INC

Spatial Autoregressive Model (SAR)

spatial_autoregressive <- lagsarlm(HOVAL ~ INC + DISCBD + CRIME +  PLUMB, data = columbus, listw = rswm_queen, Durbin = FALSE)
summary(spatial_autoregressive)
## 
## Call:lagsarlm(formula = HOVAL ~ INC + DISCBD + CRIME + PLUMB, data = columbus, 
##     listw = rswm_queen, Durbin = FALSE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.385677 -0.241482 -0.069741  0.179708  0.896593 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  2.814062   0.679476  4.1415  3.45e-05
## INC          0.108662   0.146758  0.7404 0.4590460
## DISCBD       0.092141   0.054392  1.6940 0.0902577
## CRIME       -0.012186   0.004119 -2.9585 0.0030917
## PLUMB        0.044220   0.013373  3.3067 0.0009439
## 
## Rho: 0.14571, LR test value: 0.59923, p-value: 0.43887
## Asymptotic standard error: 0.17199
##     z-value: 0.84724, p-value: 0.39686
## Wald statistic: 0.71782, p-value: 0.39686
## 
## Log likelihood: -8.855711 for lag model
## ML residual variance (sigma squared): 0.083643, (sigma: 0.28921)
## Number of observations: 49 
## Number of parameters estimated: 7 
## AIC: 31.711, (AIC for lm: 30.311)
## LM test for residual autocorrelation
## test value: 2.8762, p-value: 0.089896

Las variables muestran ser significativas por si solas (de nuevo, menos INC), pero el modelo no lo es.

Spatial Error Model (SEM)

spatial_error<-errorsarlm(HOVAL ~ INC + DISCBD + CRIME +  PLUMB, data = columbus, listw = rswm_queen, Durbin = FALSE)
summary(spatial_error)
## 
## Call:
## errorsarlm(formula = HOVAL ~ INC + DISCBD + CRIME + PLUMB, data = columbus, 
##     listw = rswm_queen, Durbin = FALSE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.415338 -0.213739 -0.082963  0.173198  0.853583 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  3.1812555  0.4659511  6.8274 8.644e-12
## INC          0.1513553  0.1429285  1.0590  0.289619
## DISCBD       0.1127691  0.0537446  2.0982  0.035884
## CRIME       -0.0125154  0.0039406 -3.1760  0.001493
## PLUMB        0.0410370  0.0135125  3.0370  0.002390
## 
## Lambda: 0.32329, LR test value: 2.2545, p-value: 0.13323
## Asymptotic standard error: 0.17559
##     z-value: 1.8412, p-value: 0.065598
## Wald statistic: 3.3899, p-value: 0.065598
## 
## Log likelihood: -8.0281 for error model
## ML residual variance (sigma squared): 0.079242, (sigma: 0.2815)
## Number of observations: 49 
## Number of parameters estimated: 7 
## AIC: 30.056, (AIC for lm: 30.311)

Al igual que el modelo anterior las variables muestran ser significativas por si solas (menos INC), pero el modelo no lo es.

Spatial Durbin Model (Modelo Global)

spatial_durbin <- lagsarlm(log(HOVAL) ~  log(INC) + DISCBD + CRIME +  log(PLUMB), data = columbus_shp, rswm_queen, type="mixed")
summary(spatial_durbin)
## 
## Call:lagsarlm(formula = log(HOVAL) ~ log(INC) + DISCBD + CRIME + log(PLUMB), 
##     data = columbus_shp, listw = rswm_queen, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.470576 -0.197194 -0.043567  0.129868  0.749509 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)     1.3849772  1.2623758  1.0971 0.272589
## log(INC)        0.4461315  0.1417657  3.1470 0.001650
## DISCBD          0.3861467  0.1183262  3.2634 0.001101
## CRIME          -0.0115876  0.0037014 -3.1306 0.001744
## log(PLUMB)      0.1849998  0.0634489  2.9157 0.003549
## lag.log(INC)    0.0412026  0.2890764  0.1425 0.886660
## lag.DISCBD     -0.2156472  0.1685096 -1.2797 0.200639
## lag.CRIME       0.0105550  0.0099407  1.0618 0.288330
## lag.log(PLUMB)  0.0546953  0.1217712  0.4492 0.653313
## 
## Rho: 0.12032, LR test value: 0.36297, p-value: 0.54686
## Asymptotic standard error: 0.19442
##     z-value: 0.61884, p-value: 0.53602
## Wald statistic: 0.38296, p-value: 0.53602
## 
## Log likelihood: -3.37305 for mixed model
## ML residual variance (sigma squared): 0.066975, (sigma: 0.2588)
## Number of observations: 49 
## Number of parameters estimated: 11 
## AIC: 28.746, (AIC for lm: 27.109)
## LM test for residual autocorrelation
## test value: 0.10366, p-value: 0.74748

Nuevamente este modelo no es significativo, sus variables sí, pero las lag de las variables, que muestra la relación con los vecinos no lo son, a pesar hasta ahora es el modelo con menor AIC.

Geographic Weighted Regression (GWR) (Modelo Local)

bw1 <- bw.gwr(log(HOVAL) ~  log(INC) + DISCBD + CRIME +  log(PLUMB), 
              approach = "AIC", adaptive = T, data = columbus_shp)
## Adaptive bandwidth (number of nearest neighbours): 37 AICc value: 36.44385 
## Adaptive bandwidth (number of nearest neighbours): 31 AICc value: 41.76016 
## Adaptive bandwidth (number of nearest neighbours): 42 AICc value: 33.03969 
## Adaptive bandwidth (number of nearest neighbours): 44 AICc value: 32.26819 
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 31.73515 
## Adaptive bandwidth (number of nearest neighbours): 47 AICc value: 31.62451 
## Adaptive bandwidth (number of nearest neighbours): 48 AICc value: 31.38235 
## Adaptive bandwidth (number of nearest neighbours): 48 AICc value: 31.38235
m.gwr <- gwr.basic(log(HOVAL) ~  log(INC) + DISCBD + CRIME +  log(PLUMB), adaptive = T, data = columbus_shp, bw = bw1) 
m.gwr
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2023-05-12 22:49:53 
##    Call:
##    gwr.basic(formula = log(HOVAL) ~ log(INC) + DISCBD + CRIME + 
##     log(PLUMB), data = columbus_shp, bw = bw1, adaptive = T)
## 
##    Dependent (y) variable:  HOVAL
##    Independent variables:  INC DISCBD CRIME PLUMB
##    Number of data points: 49
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49476 -0.25719 -0.02574  0.19306  0.91011 
## 
##    Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)  2.50278    0.53624   4.667 2.87e-05 ***
##    log(INC)     0.39832    0.15766   2.526 0.015197 *  
##    DISCBD       0.15130    0.05474   2.764 0.008311 ** 
##    CRIME       -0.01200    0.00428  -2.804 0.007485 ** 
##    log(PLUMB)   0.22453    0.06183   3.631 0.000732 ***
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 0.3005 on 44 degrees of freedom
##    Multiple R-squared: 0.5572
##    Adjusted R-squared: 0.517 
##    F-statistic: 13.84 on 4 and 44 DF,  p-value: 2.18e-07 
##    ***Extra Diagnostic information
##    Residual sum of squares: 3.974458
##    Sigma(hat): 0.2907971
##    AIC:  27.97132
##    AICc:  29.97132
##    BIC:  13.67316
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: bisquare 
##    Adaptive bandwidth: 48 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Euclidean distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                    Min.   1st Qu.    Median   3rd Qu.    Max.
##    Intercept   1.755903  2.162521  2.366178  2.434580  2.6256
##    log(INC)    0.349911  0.396317  0.444324  0.515947  0.6390
##    DISCBD      0.105533  0.135897  0.152644  0.187936  0.2366
##    CRIME      -0.014344 -0.013254 -0.011869 -0.010291 -0.0059
##    log(PLUMB)  0.157767  0.205977  0.234621  0.262730  0.3086
##    ************************Diagnostic information*************************
##    Number of data points: 49 
##    Effective number of parameters (2trace(S) - trace(S'S)): 10.15603 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 38.84397 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 31.38235 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 15.84339 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): -8.766678 
##    Residual sum of squares: 3.337269 
##    R-square value:  0.6282141 
##    Adjusted R-square value:  0.5284394 
## 
##    ***********************************************************************
##    Program stops at: 2023-05-12 22:49:53

Este muestra ser el mejor modelo al ser significativas tanto sus variables con como el modelo mismo, incluso es el modelo con menor AIC.

gwr_sf = st_as_sf(m.gwr$SDF)
gwr_sf$exp_residuals <- exp(gwr_sf$residual)
tm_shape(gwr_sf) +
  tm_polygons(col = "exp_residuals", palette="PuBu", style="quantile", n=8, title="Residuals") +
  tm_layout(title= 'Regression Residuals',  title.position = c('right', 'top'))

moran.test(gwr_sf$exp_residuals, rswm_queen) 
## 
##  Moran I test under randomisation
## 
## data:  gwr_sf$exp_residuals  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.26187, p-value = 0.3967
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.001826498      -0.020833333       0.007487432

No se muestra autocorrelación en lo residuales, pues el Moran test aparte de salir bajo salió no significativo, esto muestra que es un buen modelo.

Diagnóstico de Resultados Estimados

Multicolinealidad

vif(non_spatial_model)
##      INC   DISCBD    CRIME    PLUMB 
## 1.882200 2.835448 2.720648 1.541956

El menor valor posible de VIF es uno es decir que hay ausencia de multicolinealidad. Como regla general, un valor VIF que mayo 5 o 10 indica una cantidad problemática de colinealidad (James et al. 2014), en este caso ninguna de las variables independientes muestra multicolinearidad.

Lagrange Multiplier Diagnostic for Spatial Dependence (LMlag)

lm.LMtests(non_spatial_model,rswm_queen,test=c("RLMlag"))
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = HOVAL ~ INC + DISCBD + CRIME + PLUMB, data =
## columbus)
## weights: rswm_queen
## 
## RLMlag = 0.982, df = 1, p-value = 0.3217

Al resultar la prueba no significativa, la prueba LM de dependencia espacial nos dice que no debemos considerar el spatial lag de la variable dependiente en nuestro modelo de regresión.

Lagrange Multiplier Diagnostic for Spatial Error Dependence (LMerr)

lm.LMtests(non_spatial_model,rswm_queen,test=c("RLMerr"))
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = HOVAL ~ INC + DISCBD + CRIME + PLUMB, data =
## columbus)
## weights: rswm_queen
## 
## RLMerr = 2.4072, df = 1, p-value = 0.1208

Nuevamente con la prueba siendo no significativa lpara la dependencia del error espacial, no debemos considerar el spattial lag del término de error en nuestro modelo de regresión.

Autocorrelación Espacial de los residuales estimados (εi)

# detectar residuos de regresión espacialmente autocorrelacionados / no espacialmente autocorrelacionados
columbus_shp$non_spatial_regression_residuals <-non_spatial_model$residuals
moran.test(non_spatial_model$residuals, rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  non_spatial_model$residuals  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 1.7113, p-value = 0.04352
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.136823894      -0.020833333       0.008487864

Ya que estadística el resultado de la prueba es bajo y significativo, podemos decir que no hay concentración de residuales en Columbus. Este resultado significa que la especificación de nuestro modelo de regresión podría considerarse como una especificación correcta.

gwr_sf$exp_residuals <- exp(gwr_sf$residual)

moran.test(exp(spatial_autoregressive$residuals), rswm_queen) 
## 
##  Moran I test under randomisation
## 
## data:  exp(spatial_autoregressive$residuals)  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 1.0029, p-value = 0.158
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.068148788      -0.020833333       0.007872207
moran.test(exp(spatial_error$residuals), rswm_queen) 
## 
##  Moran I test under randomisation
## 
## data:  exp(spatial_error$residuals)  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.36551, p-value = 0.3574
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.011831447      -0.020833333       0.007986536
moran.test(exp(spatial_durbin$residuals), rswm_queen) 
## 
##  Moran I test under randomisation
## 
## data:  exp(spatial_durbin$residuals)  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.18907, p-value = 0.425
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.003796453      -0.020833333       0.008119778
moran.test(gwr_sf$exp_residuals, rswm_queen)
## 
##  Moran I test under randomisation
## 
## data:  gwr_sf$exp_residuals  
## weights: rswm_queen    
## 
## Moran I statistic standard deviate = 0.26187, p-value = 0.3967
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.001826498      -0.020833333       0.007487432

Ninguno de los modelos salió significativo por lo que podemos decir que no exitse correlación en los residuales del resto de modelos.

Selección de Modelo

Especificar e interpretar criterio de selección de modelo

export_summs(non_spatial_model, spatial_autoregressive, spatial_error, spatial_durbin)
Model 1Model 2Model 3Model 4
(Intercept)3.24 ***2.81 ***3.18 ***1.38   
(0.51)   (0.68)   (0.47)   (1.26)  
INC0.11    0.11    0.15          
(0.16)   (0.15)   (0.14)         
DISCBD0.11 *  0.09    0.11 *  0.39 **
(0.05)   (0.05)   (0.05)   (0.12)  
CRIME-0.01 *  -0.01 ** -0.01 ** -0.01 **
(0.00)   (0.00)   (0.00)   (0.00)  
PLUMB0.05 ** 0.04 ***0.04 **       
(0.01)   (0.01)   (0.01)         
rho       0.15           0.12   
       (0.17)          (0.19)  
lambda              0.32          
              (0.18)         
log(INC)                     0.45 **
                     (0.14)  
log(PLUMB)                     0.18 **
                     (0.06)  
lag.log(INC)                     0.04   
                     (0.29)  
lag.DISCBD                     -0.22   
                     (0.17)  
lag.CRIME                     0.01   
                     (0.01)  
lag.log(PLUMB)                     0.05   
                     (0.12)  
N49       49       49       49      
R20.54    0.54    0.57    0.63   
*** p < 0.001; ** p < 0.01; * p < 0.05.

Tomando en cuenta toda la información anterior de los residuos y el AIC el modelo a elegir será el cuarto, pues fue el que tuvo e menor AIC y sí tuvo significancia.

Describir los principales 5-7 hallazgos identificados a partir de los resultados de ESDA y del modelo seleccionado

  • Se muestra significacia espacial tanto para INC como PLUMB
  • Ninguno de los modelos de regresión usados muestra autocorrelación espacial en los residuales, por lo que podemos decir que en general los modelos pueden ser apropiada para predecir el valor de una casa
  • El crimen muestra ser significativo y reducir el valor de las casasa
  • El la falta de plomería de counties vecinos afecta en 0.18 el valor de las casas de un countie en particular
  • El ingreso de counties vecinos afecta en 0.45 el valor de las casas de un countie en particular

Visualizar e interpretar a través de mapa la predicción de los valores de la principal variable de interés (variable dependiente)

A continuación se mostrarán las predicciones locales de la variable dependiente y de las idenpendientes que se mostraron significativas

gwr_sf$y_predicted <- exp(gwr_sf$yhat)
tm_shape(gwr_sf) +
  tm_polygons(col = "y_predicted", palette="PuBu", style="quantile", n=8) +
   tm_layout(title= 'HOVAL',  title.position = c('right', 'top'))

Para la variable a predecir se pude ver como se predice un mayor valor de vivienda en las extremos Columbus, siendo el este y el norte las regiones con mayor valor de vivienda (HOVAL).

tm_shape(gwr_sf) +
  tm_polygons(col = "DISCBD", palette="PuBu", style="quantile", n=8, title="t-statistic") +
  tm_layout(title= 'DISCBD',  title.position = c('right', 'top'))

tm_shape(gwr_sf) +
  tm_polygons(col = "CRIME", palette="PuBu", style="quantile", n=8, title="t-statistic") +
  tm_layout(title= 'CRIME',  title.position = c('right', 'top'))

Tanto para el crimen como para la distancia se muestra concentración en la parte norte de Columbus.

Refrencias

Abdishakur. (2019, November 5). What is Exploratory Spatial Data Analysis (ESDA)? - Towards Data Science. Medium; Towards Data Science. https://towardsdatascience.com/what-is-exploratory-spatial-data-analysis-esda-335da79026ee

ArcGIS. (2019). Conceptos básicos del análisis de regresión—Ayuda | Documentación. Arcgis.com. https://desktop.arcgis.com/es/arcmap/10.7/tools/spatial-statistics-toolbox/regression-analysis-basics.htm

Bandwidth selection for basic GWR. (2022). bw.gwr function - RDocumentation. Rdocumentation.org. https://www.rdocumentation.org/packages/GWmodel/versions/2.2-9/topics/bw.gwr

Dall’erba, S. (2009). Exploratory Spatial Data Analysis. International Encyclopedia of Human Geography, 683–690. https://doi.org/10.1016/b978-008044910-4.00433-8

GIS&T. (2016). spatial non-stationarity | GIS&T Body of Knowledge. Ucgis.org. https://gistbok.ucgis.org/topic-keywords/spatial-non-stationarity

GISGeography. (2014, July 15). Spatial Autocorrelation and Moran’s I in GIS. GIS Geography. https://gisgeography.com/spatial-autocorrelation-moran-i-gis/

Guélat, J. (2013). Spatial autocorrelation (modelling). Amazonaws.com. http://rstudio-pubs-static.s3.amazonaws.com/9687_cc323b60e5d542449563ff1142163f05.html

James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. (2014). An Introduction to Statistical Learning: With Applications in R. Springer Publishing Company, Incorporated.

R: Geographically Weighted Random Forest. (2019). R-Project.org. https://search.r-project.org/CRAN/refmans/SpatialML/html/grf.html

R Spatial. (2019). Spatial autocorrelation —. Rspatial.org. https://rspatial.org/analysis/3-spauto.html